#include "Modeling\CUDADataStruct.h"
//#include ".\Modeling\RenderBaseStruct.h"
#include ".\CUDA\utility\inc\cutil.h"
#include ".\CUDA\utility\inc\cudpp\cudpp.h"
#include <stdio.h>

typedef unsigned int uint;

#define UNIFORM_KNOT 1

#define NTHREADS		128
#define PIXELTHREADS	64
//const int nThreads = 128;


#define DEBUG_BUFFERS 0
extern void dumpBuffer(uint *d_buffer, int nelements);
extern void dumpBufferfloat(float *d_buffer, int nelements);

const uint pixelBits = 10;
uint patchBits, pixelMask, patchMask;


uint3 deg, degShift, seg, segShift;
uint degSize, segSize;
float3 boxSize, box0;

uint nPatch;

#ifndef UNIFORM_KNOT
__constant__ float* dBoxX = 0;
__constant__ float* dBoxY = 0;
__constant__ float* dBoxZ = 0;
#endif

__constant__ int* dPatchIndex = 0;
__constant__ float* dFuncData = 0;
__device__	 float* dPolarData = 0;
texture<float, 1, cudaReadModeElementType> polarTex;


//the variable below are defined in the way described in "Progress 06.22"
uint maxSeg;
uint maxSilPoint;
uint actualSilPoint;
//uint maxPolarIndexBlock;
//uint nPatchEveryBlock;


size_t* dMaxSilPoint  = NULL;
__device__ uint* dSilOccupied = 0;
__device__ uint* dSilOccupiedCompact = 0;
__device__ uint* dSilOccupiedValid = 0;
__device__ uint* dSilOccupiedSort = 0;
__device__ uint* dSilOccupiedSortPad = 0;
//__device__ uint* dSilOccupiedSortIndex = 0;

CUDPPHandle plan;

__constant__ MULTI_SIL_RENDER* dMultiSilRender;
int renderWidth, renderHeight;
//float4 viewProj0, viewProj1, viewProj2, viewProj3;
float4 invViewProj0, invViewProj1, invViewProj2, invViewProj3;
//float3 eye;
//float errorXTol;
//LIGHT light[3];
//MATERIAL material;
//int nLight;

__device__ void normalize(float* a)
{
	float sum;
	sum=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
	if(sum > 1.0e-10f)
    {
		a[0]/=sum; a[1]/=sum; a[2]/=sum;
	}
}

__device__ void addambient(MATERIAL *mat, float rgb[])
{
	rgb[0]+=mat->ka*mat->d[0]*mat->a[0];
	rgb[1]+=mat->ka*mat->d[1]*mat->a[1];
	rgb[2]+=mat->ka*mat->d[2]*mat->a[2];
}

__device__ void addPhong(MATERIAL *mat, LIGHT* light, float norm[], float ray[], float rgb[])
{
	float H[3] ={(-light->dir[0]-ray[0])/2.0f, (-light->dir[1]-ray[1])/2.0f, (-light->dir[2]-ray[2])/2.0f};
	normalize(H);
	float phong = norm[0]*H[0]+norm[1]*H[1]+norm[2]*H[2];
	phong *= phong;
    phong *= phong;
    phong *= mat->ks;
	rgb[0]+=phong*light->rgb[0];
	rgb[1]+=phong*light->rgb[1];
	rgb[2]+=phong*light->rgb[2];
}

__device__ void addDiffuse(MATERIAL *mat, LIGHT* light, float diff, float norm[], float rgb[])
{
	rgb[0]+=mat->kd*diff*light->rgb[0];
	rgb[1]+=mat->kd*diff*light->rgb[1];
	rgb[2]+=mat->kd*diff*light->rgb[2];
}

__device__ void GaussEliminationGPU(float A[][3], float b[])
{
	int i,j,k,ik;
	float mik, temp;
	for(k=0; k<3; k++)
	{
		mik=-1.0f;
		for(i=k; i<3; i++)
		{
			if(fabs(A[i][k])>mik)
			{
				mik=fabs(A[i][k]);
				ik=i;
			}
		}
		if (mik<1.0e-10f)
		{
			b[0]=b[1]=b[2]=0.0f;
			return;
		}
		for(j=k; j<3; j++)
		{
			//swap(A[ik][j], A[k][j]);
			temp = A[ik][j];
			A[ik][j]=A[k][j];
			A[k][j]=temp;
		}
		//swap(b[ik], b[k]);
		temp = b[ik];
		b[ik] = b[k];
		b[k] = temp;
		b[k] /= A[k][k];
		for(i=3-1; i>=k; i--)
			A[k][i]/=A[k][k];
		for(i=k+1; i<3; i++)
		{
			b[i] -= A[i][k]*b[k];
			for (j=3-1; j>k; j--)
				A[i][j] -= A[i][k]*A[k][j];
		}
	}
	for(i=3-1; i>=0; i--)
	{
		for (j=i+1; j<3; j++)
			b[i] -= A[i][j]*b[j];
	}
}

__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}

__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossomPolarTex(float u, float v, float w, int nx, int ny, int nz, int index, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += tex1Dfetch(polarTex, index++) * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossomNormal(float u, float v, float w, int nx, int ny, int nz, uint3 degShift, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}
	func += degShift.x;
	GetBernstein(nn, ny, v);
	for (int i=0; i<nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx-1, u);
	re[0] = 0.0f;
	for (int i=0; i<nx; i++)
	{
		re[0] += f2[i]*nn[i];
	}

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<=nz; k++)
				f1[i][j] += *func++ * nn[k];
		}
		func += degShift.y;
	}
	GetBernstein(nn, ny-1, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[1] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[1] += f2[i]*nn[i];
	}

	//Z direction:
	GetBernstein(nn, nz-1, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<=ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<nz; k++)
				f1[i][j] += *func++ * nn[k];
			func++;
		}
	}
	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[2] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[2] += f2[i]*nn[i];
	}
}



extern "C"
void InitMulitSilSurfaceCUDA(float* pFuncData, int* Index, int* indexInverse, int NumPatch, MULTI_SIL_SURFACE* multiSilSurface, int surfWidth, int surfHeight)
{
	patchBits = 8*sizeof(uint)-2*pixelBits;
	uint maxNumPatch = 1<<patchBits;
	if (NumPatch>maxNumPatch)
	{
		printf("Too many patchs!\n"); 
		return;
	}

	pixelMask = (1<<pixelBits)-1;
	patchMask = (1<<patchBits)-1;


	nPatch = NumPatch;
	deg.x=multiSilSurface->degX; deg.y=multiSilSurface->degY; deg.z=multiSilSurface->degZ;
	seg.x=multiSilSurface->segX; seg.y=multiSilSurface->segY; seg.z=multiSilSurface->segZ;
	degShift.x = (deg.y+1)*(deg.z+1); degShift.y = deg.z+1; degShift.z = 1;
	segShift.x = seg.y*seg.z; segShift.y = seg.z; segShift.z = 1;
	maxSeg = max(seg.x, seg.y); maxSeg = max(maxSeg, seg.z);
	degSize = (deg.x+1)*degShift.x; segSize = seg.x*segShift.x;

	int funcSize = NumPatch*degSize*sizeof(float);	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dFuncData, funcSize) );
	CUDA_SAFE_CALL( cudaMemcpy(dFuncData, pFuncData, funcSize, cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dFuncData, pFuncData, size);	

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPatchIndex, segSize*sizeof(int)) );
	CUDA_SAFE_CALL( cudaMemcpy(dPatchIndex, Index, segSize*sizeof(int), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dPatchIndex, Index, segX*segY*segZ*sizeof(uint));

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPolarData, nPatch*degSize*sizeof(float)) );
	//cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
	CUDA_SAFE_CALL( cudaBindTexture(0, polarTex, dPolarData, cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat)/*channelDesc*/) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dMultiSilRender, sizeof(MULTI_SIL_RENDER)) );
	
	renderWidth = surfWidth;
	renderHeight = surfHeight;
 	int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupied, silSize) );		
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedValid, silSize) );
	maxSilPoint = 3*renderWidth*maxSeg;
	//maxSilPoint = renderHeight*renderWidth;
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedCompact, maxSilPoint*sizeof(uint)) );
    CUDA_SAFE_CALL( cudaMalloc( (void**) &dMaxSilPoint, sizeof(size_t)));

#if	UNIFORM_KNOT	
	box0.x = multiSilSurface->boxX[0];
	box0.y = multiSilSurface->boxY[0];
	box0.z = multiSilSurface->boxZ[0];
	boxSize.x = multiSilSurface->boxX[1]-multiSilSurface->boxX[0];
	boxSize.y = multiSilSurface->boxY[1]-multiSilSurface->boxY[0];
	boxSize.z = multiSilSurface->boxZ[1]-multiSilSurface->boxZ[0];
#else	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxX, (seg.x+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxY, (seg.y+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxZ, (seg.z+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float), cudaMemcpyHostToDevice) );	
	CUDA_SAFE_CALL( cudaMemcpy(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float), cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float));
#endif
}


__global__ void GetPolarData(float* FuncData, float* PolarData, int* PatchIndex, float3 eye, uint3 deg, uint3 degShift, uint3 segShift, uint funcSize, float3 box0, float3 boxSize)
{
	int patchIndex = PatchIndex[blockIdx.x];
	if(patchIndex<0)		return;


	extern __shared__ float funcPolar[];
	float* func = funcPolar;
	float* polar = &func[funcSize];
	float* temp = &polar[funcSize];

	//__shared__ float func[27];
	//__shared__ float polar[27];
	//__shared__ float temp[27];

	//Degree elevation:
	float A[3][2] = {{1.0f,0.0f}, {0.5f, 0.5f}, {0.0f,1.0f}};

	uint ind = threadIdx.x;	
	uint polarIndex = blockIdx.x;
	func[ind] = FuncData[patchIndex*funcSize+ind];
	polar[ind] = 0.0f;
	temp[ind] = 0.0f;
	__syncthreads();
		
	/**/
	uint patchIndexX = polarIndex/segShift.x;
	float patchX = patchIndexX*boxSize.x + box0.x;
	float eyeU = (eye.x-patchX)/boxSize.x;
	uint indexX = ind/degShift.x;
	if(indexX<deg.x)
		temp[ind] = deg.x*(eyeU*func[ind+degShift.x] + (1-eyeU)*func[ind]);
	__syncthreads();
	polar[ind]  += A[indexX][0]*temp[ind-indexX*degShift.x] + A[indexX][1]*temp[ind-indexX*degShift.x + degShift.x]; //+A[indexX][2]*polarJacobi[...+2*degShift.x]
	

	
	
	uint patchIndexY = polarIndex%segShift.x/segShift.y;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float eyeV = (eye.y-patchY)/boxSize.y;
	uint indexY = ind%degShift.x/degShift.y;
	if(indexY<deg.y)
		temp[ind] = deg.y*(eyeV*func[ind+degShift.y] + (1-eyeV)*func[ind]);
	__syncthreads();
	polar[ind] += A[indexY][0]*temp[ind-indexY*degShift.y] + A[indexY][1]*temp[ind-indexY*degShift.y + degShift.y]; //+A[indexY][2]*polarJacobi[...+2*degShift.y]
	
	uint patchIndexZ = polarIndex%segShift.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;
	float eyeW = (eye.z-patchZ)/boxSize.z;
	uint indexZ = ind%degShift.y;
	if(indexZ<deg.z)
		temp[ind] = deg.z*(eyeW*func[ind+degShift.z] + (1-eyeW)*func[ind]);
	__syncthreads();
	polar[ind] += A[indexZ][0]*temp[ind-indexZ*degShift.z] + A[indexZ][1]*temp[ind-indexZ*degShift.z + degShift.z]; //+A[indexZ][2]*polarJacobi[...+2*degShift.z]
	/**/

	PolarData[patchIndex*funcSize+ind] = polar[ind];	
}

__global__ void NoSortCheckSilPointUniform(unsigned char* surface, size_t pitch, int width, int height, uint pixelBits, uint* dSilOccupied, uint* dSilOccupiedValid)
{	
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	int y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=width || y>=height)		return;
	uint pixelI = y*width+x;
    float* pixel = (float*)(surface + y*pitch) + 4*x;
	if( pixel[0] == 0.0f || pixel[0] == 1.0f)		return;

	dSilOccupied[pixelI] = ((x<<pixelBits)|y);
	dSilOccupiedValid[pixelI] = 1;
}

__global__ void NoSortUpdateSilPointKernel(unsigned char*surfaceIn, unsigned char*surfaceOut, size_t pitch,
											MULTI_SIL_RENDER* dMultiSilRender, uint pixelBits,
											uint* SilOccupiedCompact, float* FuncData, int* PatchIndex,
											uint3 deg, uint degSize, uint3 seg, uint3 segShift,
											float3 boxSize, float3 box0//,float* elimation
										   )
{
	uint pixelX = SilOccupiedCompact[threadIdx.x + blockIdx.x*blockDim.x];
	uint pixelY = pixelX<<(sizeof(uint)*8-pixelBits)>>(sizeof(uint)*8-pixelBits);
	pixelX = pixelX>>pixelBits;
	float* pixelInfo = ((float*)(surfaceIn+pixelY*pitch) + pixelX*4 +2);

	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *(pixelInfo);

	float xyz[3] = {0.0f, 0.0f, 0.0f};

	float worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	int patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	int patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	int patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;

	int patchIndex;
	float errF = 0.0f;
	float deltaPoint = 0.0f;
	float delta = 0.01f;
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	int stepX=0, stepY=0, stepZ=0;
	
	int pixelRefine, pixelIter;
	int dir; 
	float pad;
	float* pixel;
	
	float plane[4] = {0.0f, 0.0f, 0.0f, 0.0f};

	//if(pixelInfo.x<0.3f)
	{
		plane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
		plane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
		plane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
		plane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	}
	//else
	//{
	//	plane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	//	plane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	//	plane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	//	plane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	//}
	plane[0] *= boxSize.x;
	plane[1] *= boxSize.y;
	plane[2] *= boxSize.z;

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];
		
		stepX = floorf(xyz[0]);
		stepY = floorf(xyz[1]);
		stepZ = floorf(xyz[2]);

		patchIndexX += stepX;
		patchIndexY += stepY;
		patchIndexZ += stepZ;

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		patchIndex *= degSize;
		
		patchX += stepX*boxSize.x;
		patchY += stepY*boxSize.y;
		patchZ += stepZ*boxSize.z;

		plane[3] += stepX*plane[0] + stepY*plane[1] + stepZ*plane[2];

		xyz[0] -= stepX;
		xyz[1] -= stepY;
		xyz[2] -= stepZ;

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &fValue[1]);
		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];

		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		EvalBezierBlossomPolarTex(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][1]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, patchIndex, &jacobi[1][2]);

		jacobi[1][0] = (jacobi[1][0]-fValue[1])/delta;
		jacobi[1][1] = (jacobi[1][1]-fValue[1])/delta;
		jacobi[1][2] = (jacobi[1][2]-fValue[1])/delta;

		jacobi[2][0] = plane[0];
		jacobi[2][1] = plane[1];
		jacobi[2][2] = plane[2];


		//elimation[eliIndex  ] = jacobi[0][0];
		//elimation[eliIndex+1] = jacobi[0][1];
		//elimation[eliIndex+2] = jacobi[0][2];
		//elimation[eliIndex+3] = jacobi[1][0];
		//elimation[eliIndex+4] = jacobi[1][1];
		//elimation[eliIndex+5] = jacobi[1][2];
		//elimation[eliIndex+6] = plane[0];
		//elimation[eliIndex+7] = plane[1];
		//elimation[eliIndex+8] = plane[2];
		//elimation[eliIndex+9] = fValue[0];
		//elimation[eliIndex+10] = fValue[1];
		//elimation[eliIndex+11] = fValue[2];
		//patchIndex+=15;
		//for(int i=0; i<15; i++)
		//	elimation[eliIndex+i] = tex1Dfetch(polarTex, patchIndex++);

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= dMultiSilRender->errorTol)					break;

		GaussEliminationGPU(jacobi, fValue);

		//elimation[eliIndex+12] = fValue[0];
		//elimation[eliIndex+13] = fValue[1];
		//elimation[eliIndex+14] = fValue[2];


		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= dMultiSilRender->errorXTol)			break;	
	}

	//if(pixelY==115&&pixelX==491)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000000;
	//	return;
	//}
	//else
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 0;
	//	return;
	//}

	//if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 0;
	//	return;
	//}

	//if(patchIndex<0)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1;
	//	return;
	//}
	//	
	//if(errF>dMultiSilRender->errorTol/* && deltaPoint>dMultiSilRender->errorXTol*/)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 2;
	//	return;
	//}

	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*100 && deltaPoint<dMultiSilRender->errorXTol) )
	{
		xyz[0] = xyz[0]*boxSize.x+patchX;
		xyz[1] = xyz[1]*boxSize.y+patchY;
		xyz[2] = xyz[2]*boxSize.z+patchZ;

		fValue[0] = xyz[0]*dMultiSilRender->viewProj[0][0] + xyz[1]*dMultiSilRender->viewProj[1][0] + xyz[2]*dMultiSilRender->viewProj[2][0] + dMultiSilRender->viewProj[3][0]; 
		//fValue[1] = xyz[0]*dMultiSilRender->viewProj[0][1] + xyz[1]*dMultiSilRender->viewProj[1][1] + xyz[2]*dMultiSilRender->viewProj[2][1] + dMultiSilRender->viewProj[3][1]; 
		fValue[2] = xyz[0]*dMultiSilRender->viewProj[0][3] + xyz[1]*dMultiSilRender->viewProj[1][3] + xyz[2]*dMultiSilRender->viewProj[2][3] + dMultiSilRender->viewProj[3][3]; 


		//if(pixelInfo.x<0.3f)
		{
			pixelRefine = floorf(((1.0f+fValue[0]/fValue[2])*dMultiSilRender->renderWidth)/2.0f);
			pixelRefine = min(pixelRefine, dMultiSilRender->renderWidth-1);
			pixelRefine = max(0, pixelRefine);
/*
			pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelRefine);
			(*pixel++) = 0.0f;
			(*pixel++) = 0.0f;
			(*pixel++) = 0.0f;
			int aaa = (patchIndexX+patchIndexY+patchIndexZ)%3;
			pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelRefine + aaa);
			(*pixel) = 1.0f;
*/
			pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelX + 2);
			dir = (pixelRefine>pixelX) ? 1 : -1;
			pad = *(pixelInfo-4*dir);
			//pad = 0.0f;
			for(pixelIter=pixelX; (pixelIter!=pixelRefine/*&&pixelIter>=0&&pixelIter<dMultiSilRender->renderWidth*/); pixelIter+=dir, pixel+=4*dir)
				(*pixel) = pad;

			(*pixel--) = 0.0f;
			(*pixel--) = 0.0f;
			(*pixel--) = 1.0f;
			pixel-=4*dir;
			(*pixel++) = 1.0f;
			(*pixel++) = 0.0f;
			(*pixel++) = 0.0f;

			SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000 + (pixelRefine-pixelX);
			
		}
		//else
		//{
		//	pixelRefine = floorf(((1.0f-fValue[1]/fValue[2])*dMultiSilRender->renderHeight)/2.0f);
		//	dir = (pixelRefine>pixelY) ? 1 : -1;
		//	pad = (pixelRefine>pixelY) ? pixelInfo.y : pixelInfo.w;		
		//	//pad = 0.0f;
		//	for(pixelIter=pixelY; (pixelIter!=pixelRefine&&pixelIter>=0&&pixelIter<dMultiSilRender->renderHeight); pixelIter+=dir, pixel+=(pitch/4)*dir)
		//		(*pixel) = pad;
		//	if(pixelIter==pixelRefine)
		//		(*pixel) = pad;
		//	//SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = pixelRefine-pixelY;
		//}
	}


	

	//if(errF<dMultiSilRender->errorTol || deltaPoint<dMultiSilRender->errorXTol)
	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	//{
	//	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	//	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	//	xyz[2] = (xyz[2]-patchZ)/boxSize.z;
	//}
	//else
	//{
	//	xyz[0] += fValue[0];
	//	xyz[1] += fValue[0];
	//	xyz[2] += fValue[0];
	//	patchIndexX -= stepX;
	//	patchIndexY -= stepY;
	//	patchIndexZ -= stepZ;
	//}
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	patchX = patchIndexX*boxSize.x + box0.x;
	patchY = patchIndexY*boxSize.y + box0.y;
	patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;
	

	plane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	plane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	plane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	plane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	plane[0] *= boxSize.x;
	plane[1] *= boxSize.y;
	plane[2] *= boxSize.z;
	fValue[0] = fValue[1] = fValue[2] = 0.0f;

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];
		
		stepX = floorf(xyz[0]);
		stepY = floorf(xyz[1]);
		stepZ = floorf(xyz[2]);

		patchIndexX += stepX;
		patchIndexY += stepY;
		patchIndexZ += stepZ;

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		patchIndex *= degSize;
		
		patchX += stepX*boxSize.x;
		patchY += stepY*boxSize.y;
		patchZ += stepZ*boxSize.z;

		plane[3] += stepX*plane[0] + stepY*plane[1] + stepZ*plane[2];

		xyz[0] -= stepX;
		xyz[1] -= stepY;
		xyz[2] -= stepZ;

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &fValue[1]);
		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];


		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		EvalBezierBlossomPolarTex(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][1]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, patchIndex, &jacobi[1][2]);

		jacobi[1][0] = (jacobi[1][0]-fValue[1])/delta;
		jacobi[1][1] = (jacobi[1][1]-fValue[1])/delta;
		jacobi[1][2] = (jacobi[1][2]-fValue[1])/delta;

		jacobi[2][0] = plane[0];
		jacobi[2][1] = plane[1];
		jacobi[2][2] = plane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= dMultiSilRender->errorTol)					break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= dMultiSilRender->errorXTol)			break;	
	}
	


	//isOK=true;
	//if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 90;
	//	isOK=false;
	//	//return;
	//}

	//if(patchIndex<0)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 4;
	//	isOK=false;
	//	//return;
	//}
	//	
	//if(errF>dMultiSilRender->errorTol && deltaPoint>dMultiSilRender->errorXTol)
	//{
	//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 100;
	//	isOK=false;
	//	//return;
	//}

	

	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	//if(errF<dMultiSilRender->errorTol)
	if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*100 && deltaPoint<dMultiSilRender->errorXTol) )
	{
		xyz[0] = xyz[0]*boxSize.x+patchX;
		xyz[1] = xyz[1]*boxSize.y+patchY;
		xyz[2] = xyz[2]*boxSize.z+patchZ;

		//fValue[0] = xyz[0]*dMultiSilRender->viewProj[0][0] + xyz[1]*dMultiSilRender->viewProj[1][0] + xyz[2]*dMultiSilRender->viewProj[2][0] + dMultiSilRender->viewProj[3][0]; 
		fValue[1] = xyz[0]*dMultiSilRender->viewProj[0][1] + xyz[1]*dMultiSilRender->viewProj[1][1] + xyz[2]*dMultiSilRender->viewProj[2][1] + dMultiSilRender->viewProj[3][1]; 
		fValue[2] = xyz[0]*dMultiSilRender->viewProj[0][3] + xyz[1]*dMultiSilRender->viewProj[1][3] + xyz[2]*dMultiSilRender->viewProj[2][3] + dMultiSilRender->viewProj[3][3]; 

		//if(pixelInfo.x<0.3f)
		//{
		//	pixelRefine = floorf(((1.0f+fValue[0]/fValue[2])*dMultiSilRender->renderWidth)/2.0f);
		//	dir = (pixelRefine>pixelX) ? 1 : -1;
		//	pad = (pixelRefine>pixelX) ? *(pixelInfo-2) : *(pixelInfo+6);
		//	//pad = 1.0f;
		//	for(pixelIter=pixelX; (pixelIter!=pixelRefine&&pixelIter>=0&&pixelIter<dMultiSilRender->renderWidth); pixelIter+=dir, pixel+=4*dir)
		//		(*pixel) = pad;
		//	if(pixelIter==pixelRefine)
		//		(*pixel) = pad;
		//	SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000 + (pixelRefine-pixelX);
		//}
		//else
		{
			pixelRefine = floorf(((1.0f-fValue[1]/fValue[2])*dMultiSilRender->renderHeight)/2.0f);
			pixelRefine = min(pixelRefine, dMultiSilRender->renderHeight-1);
			pixelRefine = max(0, pixelRefine);
			
			pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelX+2);
			dir = (pixelRefine>pixelY) ? 1 : -1;
			pad = *(pixelInfo-pitch/4*dir);		
			//pad = 0.0f;
			for(pixelIter=pixelY; (pixelIter!=pixelRefine/*&&pixelIter>=0&&pixelIter<dMultiSilRender->renderHeight*/); pixelIter+=dir, pixel+=(pitch/4)*dir)
				(*pixel) = pad;

			(*pixel--) = 0.0f;
			(*pixel--) = 0.0f;
			(*pixel--) = 1.0f;
			pixel-=(pitch/4)*dir;
			(*pixel++) = 1.0f;
			(*pixel++) = 0.0f;
			(*pixel++) = 0.0f;

			SilOccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 10*(pixelRefine-pixelY)-1;
		}
	}
/**/	
}

__global__ void NoSortShadingKernel(//unsigned char*surfaceIn, 
									unsigned char*surfaceOut, size_t pitch,
									int renderWidth, int renderHeight,
									MULTI_SIL_RENDER* dMultiSilRender, uint pixelBits,
									float* FuncData, int* PatchIndex,
									uint3 deg, uint degSize, uint3 seg, uint3 segShift,
									float3 boxSize, float3 box0)
{
	int pixelX = blockIdx.x*blockDim.x + threadIdx.x;
	int pixelY = blockIdx.y*blockDim.y + threadIdx.y;
	if(pixelX>=renderWidth || pixelY>=renderHeight)		return;

	float* pixelInfo = ((float*)(surfaceOut+pixelY*pitch) + pixelX*4);
	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *(pixelInfo+2);
	if(screenZ <= 0.0001f )									
	{
		//pixelInfo[0] = 0;
		//pixelInfo[1] = 0;
		//pixelInfo[2] = 0;
		return;
	}


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float horPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float verPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

	float worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	int patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	int patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	int patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;

	horPlane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
	horPlane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
	horPlane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
	horPlane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*horPlane[0] + patchY*horPlane[1] + patchZ*horPlane[2];
	horPlane[0] *= boxSize.x;
	horPlane[1] *= boxSize.y;
	horPlane[2] *= boxSize.z;
	verPlane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	verPlane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	verPlane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	verPlane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*verPlane[0] + patchY*verPlane[1] + patchZ*verPlane[2];
	verPlane[0] *= boxSize.x;
	verPlane[1] *= boxSize.y;
	verPlane[2] *= boxSize.z;

	int patchIndex;
	float errF = 10.0f;
	float deltaPoint = 0.0f;
	float delta = 0.01f;
	int stepX=0, stepY=0, stepZ=0;

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		stepX = floorf(xyz[0]);
		stepY = floorf(xyz[1]);
		stepZ = floorf(xyz[2]);

		patchIndexX += stepX;
		patchIndexY += stepY;
		patchIndexZ += stepZ;

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		patchIndex *= degSize;
		
		patchX += stepX*boxSize.x;
		patchY += stepY*boxSize.y;
		patchZ += stepZ*boxSize.z;

		horPlane[3] += stepX*horPlane[0] + stepY*horPlane[1] + stepZ*horPlane[2];
		verPlane[3] += stepX*verPlane[0] + stepY*verPlane[1] + stepZ*verPlane[2];

		xyz[0] -= stepX;
		xyz[1] -= stepY;
		xyz[2] -= stepZ;

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		fValue[1] = horPlane[0]*xyz[0] + horPlane[1]*xyz[1] + horPlane[2]*xyz[2] + horPlane[3];
		fValue[2] = verPlane[0]*xyz[0] + verPlane[1]*xyz[1] + verPlane[2]*xyz[2] + verPlane[3];


		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];

		jacobi[1][0] = horPlane[0];
		jacobi[1][1] = horPlane[1];
		jacobi[1][2] = horPlane[2];
		jacobi[2][0] = verPlane[0];
		jacobi[2][1] = verPlane[1];
		jacobi[2][2] = verPlane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= dMultiSilRender->errorTol)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabsf(fValue[0])+fabsf(fValue[1])+fabsf(fValue[2]);
		if (deltaPoint <= dMultiSilRender->errorXTol)			break;	
	}

	//if(errF<=dMultiSilRender->errorTol*10)
	if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*100 && deltaPoint<dMultiSilRender->errorXTol) )
	{
	//MATERIAL material;
	//material.ka = dMultiSilRender->material.ka;
	//material.kd = dMultiSilRender->material.kd;
	//material.ks = dMultiSilRender->material.ks;
	//material.a[0] = dMultiSilRender->material.a[0];
	//material.a[1] = dMultiSilRender->material.a[1];
	//material.a[2] = dMultiSilRender->material.a[2];
	//material.s[0] = dMultiSilRender->material.s[0];
	//material.s[1] = dMultiSilRender->material.s[1];
	//material.s[2] = dMultiSilRender->material.s[2];

	//switch ((patchIndexX+patchIndexY+patchIndexZ)%4) 
	//{
	//	case  0:
	//		material.d[0]=1.0; material.d[1]=material.d[2]=0.0;
	//		break;
	//	case  1:
	//		material.d[1]=1.0; material.d[0]=material.d[2]=0.0;
	//		break;
	//	case  3:
	//		material.d[2]=1.0; material.d[1]=material.d[0]=0.0;
	//		break;
	//	case  2:
	//		material.d[1]=0.0; material.d[2]=material.d[0]=1.0;
	//		break;
	//}

    float rgb[3]={0.0f, 0.0f, 0.0f};

	norm[0] /= boxSize.x;
	norm[1] /= boxSize.y;
	norm[2] /= boxSize.z;
	normalize(norm);
	xyz[0] = xyz[0]*boxSize.x+patchX;
	xyz[1] = xyz[1]*boxSize.y+patchY;
	xyz[2] = xyz[2]*boxSize.z+patchZ;
    

	float ray[3]={xyz[0]-dMultiSilRender->eye[0], xyz[1]-dMultiSilRender->eye[1], xyz[2]-dMultiSilRender->eye[2]};
	normalize(ray);

    float diff;
	for(int i=0; i<dMultiSilRender->nLight; i++)  
	{
        diff = -(norm[0]*dMultiSilRender->light[i].dir[0]+norm[1]*dMultiSilRender->light[i].dir[1]+norm[2]*dMultiSilRender->light[i].dir[2]);
        if (diff>0)
        {
            //addPhong(&dMultiSilRender->material, &dMultiSilRender->light[i], norm, ray, rgb);
            addDiffuse(&dMultiSilRender->material, &dMultiSilRender->light[i], diff, norm, rgb);
			//addDiffuse(&material, &dMultiSilRender->light[i], diff, norm, rgb);
        }
	}
	addambient(&dMultiSilRender->material, rgb);
	//addambient(&material, rgb);

	//pixelInfo = ((float*)(surfaceOut+pixelY*pitch) + pixelX*4);
    pixelInfo[0] = rgb[0];
    pixelInfo[1] = rgb[1];
    pixelInfo[2] = rgb[2];
	}
	else
	{
		pixelInfo[0] = 0;
		pixelInfo[1] = 0;
		pixelInfo[2] = 0;

	}

}


extern "C"
void MultiSilRender(MULTI_SIL_RENDER* multiSilRender, void* surfaceIn, void* surfaceOut, size_t pitch, int iter)
{
	CUDA_SAFE_CALL( cudaMemcpy(dMultiSilRender, multiSilRender, sizeof(MULTI_SIL_RENDER), cudaMemcpyHostToDevice) );

	float3 eye = make_float3(multiSilRender->eye[0], multiSilRender->eye[1], multiSilRender->eye[2]);
	size_t Ns = degSize*20*sizeof(float);
#if	UNIFORM_KNOT
	GetPolarData<<<segSize, degSize, Ns>>>(dFuncData, dPolarData, dPatchIndex, eye, deg, degShift, segShift, degSize, box0, boxSize);
#else
#endif


#if DEBUG_BUFFERS
	printf("dPolarData:\n");
	dumpBufferfloat(dPolarData, degSize*nPatch);
#endif

	int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMemset(dSilOccupied, 0, silSize) );
	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedValid, 0, silSize) );
	dim3 Db = dim3(16,16);
	dim3 Dg = dim3( renderWidth/Db.x+1, renderHeight/Db.y+1 );

#if	UNIFORM_KNOT	
	NoSortCheckSilPointUniform<<<Dg,Db>>>((unsigned char*)surfaceIn, pitch, renderWidth, renderHeight, 
									pixelBits, dSilOccupied, dSilOccupiedValid);	
#else
#endif

#if DEBUG_BUFFERS
	printf("dSilOccupied:\n");
	dumpBuffer(dSilOccupiedValid, renderWidth*renderHeight);
#endif


	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedCompact, 0, maxSilPoint*sizeof(uint)) );


	CUDPPConfiguration config;
	CUDPPResult result = CUDPP_SUCCESS;	

    config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, renderWidth*renderHeight, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, dSilOccupiedCompact, dMaxSilPoint, dSilOccupied, dSilOccupiedValid, renderWidth*renderHeight);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }

	actualSilPoint=0;
	CUDA_SAFE_CALL( cudaMemcpy(&actualSilPoint, dMaxSilPoint, sizeof(size_t), cudaMemcpyDeviceToHost) );

#if DEBUG_BUFFERS
	float* elimation;
	CUDA_SAFE_CALL( cudaMalloc((void**)&elimation, actualSilPoint*15*sizeof(float)) );
	printf("dSilOccupiedCompact:\n");
	dumpBuffer(dSilOccupiedCompact, actualSilPoint);
#endif

#if	UNIFORM_KNOT	
	NoSortUpdateSilPointKernel<<<actualSilPoint/PIXELTHREADS+1, PIXELTHREADS>>>
								((unsigned char*)surfaceIn, (unsigned char*)surfaceOut, pitch,
								dMultiSilRender, pixelBits,
								dSilOccupiedCompact, dFuncData, dPatchIndex,
								deg, degSize, seg, segShift,
								boxSize, box0//,elimation								
								);
#else
#endif

#if DEBUG_BUFFERS
	printf("dSilOccupiedCompact:\n");
	dumpBufferfloat(elimation, actualSilPoint*15);
	CUDA_SAFE_CALL( cudaFree(elimation) );
	dumpBuffer(dSilOccupiedCompact, actualSilPoint);
#endif



	//dim3 Db, Dg;
	Db.x = Db.y = 8;
	Dg.x = renderWidth/Db.x+1;
	Dg.y = renderHeight/Db.y+1;

#if	UNIFORM_KNOT	
	NoSortShadingKernel<<<Dg, Db>>>(//(unsigned char*)surfaceIn,
									(unsigned char*)surfaceOut, pitch,
									renderWidth, renderHeight,
									dMultiSilRender, pixelBits,
									dFuncData, dPatchIndex,
									deg, degSize, seg, segShift,
									boxSize, box0);
#else
#endif 

/**/
}

extern "C"
void CleanUpMultiRenderCUDA()
{
	CUDA_SAFE_CALL( cudaFree(dFuncData) );
	CUDA_SAFE_CALL( cudaFree(dPatchIndex) );

	CUDA_SAFE_CALL( cudaUnbindTexture(polarTex) );
	CUDA_SAFE_CALL( cudaFree(dPolarData) );
	CUDA_SAFE_CALL( cudaFree(dMultiSilRender) );

	CUDA_SAFE_CALL( cudaFree(dSilOccupiedValid) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupied) );
	CUDA_SAFE_CALL( cudaFree(dMaxSilPoint) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupiedCompact) );


#ifndef UNIFORM_KNOT	
	CUDA_SAFE_CALL( cudaFree(dBoxX) );
	CUDA_SAFE_CALL( cudaFree(dBoxY) );
	CUDA_SAFE_CALL( cudaFree(dBoxZ) );
#endif
}
